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[i]  Prior  to  2003,  there  are  two  known  cases  where  ultrarelativistic  (>  10  MeV)  electrons 
appeared  in  the  Earth’s  inner  zone  radiation  belts  in  association  with  high  speed 
interplanetary  shocks:  the  24  March  1991  and  the  less  well  studied  21  February  1994 
storms.  During  the  March  1991  event  electrons  were  injected  well  into  the  inner  zone  on  a 
timescale  of  minutes,  producing  a  new  stably  trapped  radiation  belt  population  that 
persisted  for  ~  10  years.  More  recently,  at  the  end  of  solar  cycle  23,  a  number  of  violent 
geomagnetic  disturbances  resulted  in  large  variations  in  ultrarelativistic  electrons  in  the 
inner  zone,  indicating  that  these  events  are  less  rare  than  previously  thought.  Here  we 
present  results  from  a  numerical  study  of  shock-induced  transport  and  energization  of 
outer  zone  electrons  in  the  1-7  MeV  range,  resulting  in  a  newly  formed  10-20  MeV 
electron  belt  near  L  ~  3.  Test  particle  trajectories  are  followed  in  time-dependent  fields 
from  an  MHD  magnetospheric  model  simulation  of  the  29  October  2003  storm  sudden 
commencement  (SSC)  driven  by  solar  wind  parameters  measured  at  ACE.  The  newly 
formed  belt  is  predominantly  equatorially  mirroring.  This  result  is  in  part  due  to  an  SSC 
electric  field  pulse  that  is  strongly  peaked  in  the  equatorial  plane,  preferentially 
accelerating  equatorially  mirroring  particles.  The  timescale  for  subsequent  pitch  angle 
diffusion  of  the  new  belt,  calculated  using  quasi-linear  bounce-averaged  diffusion 
coefficients,  is  in  agreement  with  the  observed  delay  in  the  appearance  of  peak  fluxes  at 
SAMPEX  in  low  Earth  orbit.  We  also  present  techniques  for  modeling  radiation  belt 
dynamics  using  test  particle  trajectories  in  MHD  fields.  Simulations  are  performed  using 
code  developed  by  the  Center  for  Integrated  Space  Weather  Modeling. 

Citation:  Kress,  B.  T.,  M.  K.  Hudson,  M.  D.  Looper,  J.  Albert,  J.  G.  Lyon,  and  C.  C.  Goodrich  (2007),  Global  MHD  test  particle 
simulations  of  >10  MeV  radiation  belt  electrons  during  storm  sudden  commencement,  J.  Geophvs.  Res.,  1/2 ,  A09215, 
doi:  1 0. 1 029/2006JA0 12218. 


1.  Introduction 

[2]  In  comparison  with  the  Earth’s  outer  zone  radiation 
belts,  sudden  large  variations  in  the  inner  zone  energetic 
particle  fluxes  ( L  £  3)  are  rare,  occurring  only  during  very 
large  geomagnetic  storms,  usually  initiated  by  coronal  mass 
ejection  (CME)  driven  interplanetary  shocks.  This  contrast 
is  especially  evident  for  very  energetic  particles  (>  10  MeV) 
which  have  long  radial  diffusion  and  loss  timescales  allow¬ 
ing  stably  trapped  electron  and  ion  belts  to  persist  for 
months  to  years.  To  illustrate,  Figure  1  shows  a  summary 
plot  of  daily  averages  of  10-20  MeV  electron  count  rates 
versus  /.-shell  (/)  from  the  Proton/E lectron  Telescope  (PET) 
on  the  Solar  Anomalous  and  Magnetospheric  Particle  Ex- 
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plorer  (SAMPEX)  spacecraft.  The  plot  shows  10  20  MeV 
electron  count  rates  from  the  start  of  the  SAMPEX  mission 
in  July  1992  through  September  2005.  The  heightened 
fluxes  near  L  ~  2,  present  at  the  beginning  of  the  SAMPEX 
mission,  are  the  remnants  of  the  24  March  1991  CME- 
driven  interplanetary  shock  that  compressed  the  magneto¬ 
pause  inside  geosynchronous  orbit  injecting  electrons  well 
into  the  inner  zone  on  a  timescale  of  minutes  [Blake  et  a/., 
1992].  The  first  significant  enhancement  following  the 
March  1991  storm  appearing  in  Figure  1  is  associated  with 
the  21  February  1994  storm.  A  strong  ground-based  SSC 
was  observed  by  low-latitude  magnetometers  for  this  event 
(K.  Shiokawa,  private  communication,  2005),  roughly  half 
the  amplitude  of  the  March  1991  SSC  which  was  the 
strongest  on  record  [Araki  et  al. ,  1997].  This  event  was 
followed  by  a  long  period  of  relatively  little  change  during 
solar  minimum,  between  solar  cycles  22  and  23,  character¬ 
ized  by  the  slow  radial  diffusion  and  loss  timescales  evident 
in  the  plot.  The  violent  geomagnetic  storms  of  October 
November  2003  mark  the  beginning  of  the  strong  activity 
characterizing  the  declining  phase  of  the  solar  cycle.  During 
the  “Halloween  storm,”  ultrarelativistic  electrons  were 
injected  well  inside  of  L  ~  3.5  producing  a  stably  trapped 


A09215 


I  of  11 


A09215 


KRESS  ET  AL.:  SIMULATIONS  OF  >10  MeV  RADIATION  BELT  ELECTRONS 


A09215 


1994  1996  1998  2000  2002  2004 

Year 


-3  0  -2.5  -2.0  -15 

log  10  electrons/sec.  10-20  MeV 

Figure  1.  Daily  averages  of  10-20  MeV  electron  count  rates  from  the  PET  instrument  on  SAMPEX 
from  July  1992  through  September  2005.  A  new  belt  injected  during  the  29  October  2003  storm  appears 
as  a  weak  enhancement  near  1  =  2  beginning  ~  February  2004. 


radiation  belt  population  that  persisted  for  several  months 
[Looper  et  al.,  2005].  This  belt  appears  in  Figure  1  as  a 
weak  enhancement  in  fluxes  near  L  =  2  beginning  ~ 
February  2004;  and  it  is  evident  in  the  expanded  Figure  4 
of  Looper  et  al.  that  weak  fluxes,  just  above  background 
level,  extend  back  in  time  to  October -November  2003. 
Looper  et  al.  attribute  the  ~4  month  delay  in  the  appearance 
of  peak  fluxes  at  SAMPEX  (in  low  Earth  orbit)  to  a  slow 
pitch  angle  diffusion  from  a  population  initially  mirroring 
near  the  equatorial  plane.  Looper  et  al.  also  note  that 
>10  MeV  electrons  injected  during  the  24  March  1991 
storm  observed  by  the  Combined  Release  and  Radiation 
Effects  Satellite  (CRRES)  (in  near  equatorial  orbit)  had  an 
equatorial  pitch  angle  distribution  strongly  peaked  near  90°, 
suggesting  that  a  similar  physical  mechanism  is  responsible 
for  the  newly  formed  belt  in  each  case.  Several  additional 
enhancements  are  seen  in  Figure  1,  after  the  beginning  of 
2004,  associated  with  November  2004,  January  2005,  and 
May  2005  storms. 

[3]  For  comparison,  Figure  2  shows  2-6  MeV  electron 
count  rates  versus  L  (bottom)  and  10-20  MeV  electron 
count  rates  versus  L  (top),  both  from  the  PET  instrument  on 
SAMPEX  (i.e.,  the  top  of  Figure  2  is  an  expanded  version  of 
Figure  1  over  the  9  month  period  following  1  October  2003). 
A  new  belt  of  2-6  MeV  electrons  appears  in  and  below  the 


SAMPEX/PET  10-20  MeV  tog1()(#/sec) 


2003.75  2004.0  2004.25  2004.5 

SAMPEX/PET  2-6  MeV  tog10(#/sec) 


Year 

Figure  2.  (top)  The  10-20  MeV  electron  count  rate  from 
the  PET  instrument  on  SAMPEX  from  1  October  2003  to 
1  July  2004.  The  vertical  red  line  at  29  October  2003  is  a 
large  solar  energetic  electron  event  [Looper  et  a/.,  2005]. 
(bottom)  The  2-6  MeV  electron  flux  from  the  PET 
instrument  on  SAMPEX  with  the  same  timescale  as  above. 
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slot  region  1-3  days  following  the  29  October  2003  storm 
onset  [Baker  et  al. ,  2004].  This  newly  formed  2-6  MeV 
belt  has  been  attributed  both  to  adiabatic  energization  and 
transport  due  to  strongly  enhanced  magnetospheric  ultra- 
low  frequency  (ULF)  waves  [Loto'aniu  et  al .,  2006]  and  to 
local  heating  by  whistler  mode  chorus  waves  [Home  et  al ., 
2005a,  2005b;  Shprits  et  al .,  2006].  The  striking  difference 
in  timescales  and  location  of  the  appearance  of  Halloween 
storm  electrons  for  these  different  energy  ranges  suggest 
that  a  different  physical  process  is  responsible  for  the 
formation  of  new  belts  in  each  case. 

[4]  A  mechanism  that  has  been  investigated  previously,  in 
connection  with  the  formation  of  new  energetic  electron  and 
ion  belts,  is  the  inward  transport  and  adiabatic  acceleration 
of  outer  zone  electrons  by  an  electric  field  pulse  launched 
when  an  interplanetary  shock  impacts  the  magnetosphere. 
The  inductive  SSC  electric  field  pulse  accompanies  an 
increase  in  the  magnetic  field  measured  by  spacecraft  and 
ground  magnetometers.  The  leading  portion  of  the  bipolar 
electric  field  pulse  can  exceed  ~100  mV/m  [Wygant  et  al ., 
1994]  and  is  predominantly  westward.  Energetic  electrons 
in  drift  resonance  with  the  azimuthally  propagating  SSC 
electric  field  pulse  can  E  x  B  drift  inward  through  several 
L-shells,  undergoing  significant  energization  in  a  fraction  of 
a  drift  period  due  to  conservation  of  the  first  adiabatic 
invariant.  Li  et  al.  [1993]  modeled  the  formation  of  a  new 
electron  belt  at  L  ~  2.5  during  the  24  March  1991  storm  by 
following  electron  guiding  centers  restricted  to  the  equato¬ 
rial  plane  in  a  pure  dipole  magnetic  field  traversed  by  an 
analytically  modeled  bipolar  electric  field  pulse.  The  model 
was  found  to  reproduce  the  observed  electron  drift  echoes 
well.  Hudson  et  al.  [1995]  used  the  same  technique  to 
model  the  formation  of  a  new  proton  radiation  belt  at  L  ~ 
2.5,  also  observed  during  the  March  1991  event  [Blake  et 
al.,  1992].  Hudson  et  al.  [1997]  used  an  MHD-guiding 
center  test  particle  simulation  of  the  March  1991  storm, 
again  restricted  to  equatorial  plane  particle  dynamics,  to 
simulate  the  formation  of  a  new  proton  belt  earthward  of 
solar  energetic  proton  penetration.  Elkington  et  al.  [2004] 
applied  similar  techniques  to  electrons,  reproducing  the 
flux  peak  in  energy  at  13  MeV  and  L  ~  2.5  which  was 
the  signature  of  prompt  injection  on  a  drift  timescale  in 
the  CRRES  measurements.  A  recent  parametric  study  by 
Gannon  et  al.  [2005]  investigated  the  shock-induced  trans¬ 
port  and  energization  of  relativistic  electrons  in  the  magne¬ 
tosphere  using  the  Li  et  al.  model.  Gannon  et  al.  studied  the 
location  and  intensity  of  new  belts  for  a  variety  of  electric 
field  pulse  propagation  velocities  and  amplitudes  in  the 
range  from  750  to  2500  km/s  and  120  to  240  mV/m  and 
found  that  a  pulse  speed  ^1200  km/s  is  required  to  produce 
significant  flux  levels  in  the  new  belt,  also  that  a  pulse 
amplitude  1 20  mV/m  is  required  to  form  a  new  belt  inside 
ofL  ~  3. 

[5]  Similarities  between  the  March  1991  event  and  sub¬ 
sequent  injections  of  >10  MeV  electrons  shown  in  Figure  1 
lead  us  to  consider  the  possibility  that  these  additional 
enhancements  are  also  prompt  shock-induced  injections 
due  to  transport  and  energization  by  the  SSC  electric  field 
pulse  associated  with  each  event.  The  21  February  1994 
event  has  been  modeled  in  a  separate  study  limited  to 
equatorial  plane  guiding  center  test  particle  dynamics 
[Hudson  et  al.,  2006],  with  conclusions  similar  to  the 


24  March  1991  event  modeled  by  Elkington  et  al.  [2002  ]. 
In  the  present  study  of  the  29  October  2003  event,  we  use  an 
MHD  magnetospheric  model  to  simulate  an  ~5  min  period 
during  and  immediately  after  the  arrival  of  the  interplane¬ 
tary  shock  that  initiated  the  storm.  Test  particle  trajectories 
are  followed  in  MHD  fields  to  investigate  the  appearance  of 
>10  MeV  electrons  at  L  ~  2.  The  inner  boundary  conditions 
used  in  the  MHD  model  at  L  ~  2,  however,  restricts  inward 
particle  transport  at/near  the  inner  boundary,  limiting  this 
study  to  regions  L  £  2.  One  question  we  address  is  why  one 
would  expect  the  resulting  new  belt  to  have  a  pitch  angle 
distribution  strongly  peaked  near  90°  as  indicated  by  obser¬ 
vations;  thus  it  is  necessary  to  simulate  three-dimensional 
(3-D)  particle  dynamics  not  restricted  to  the  equatorial 
plane. 

[6]  In  section  2  we  present  techniques  used  for  global 
modeling  of  radiation  belt  dynamics  using  three-dimensional 
test  particle  trajectories  in  MHD  fields,  including  a  novel 
method  of  obtaining  modeled  fluxes  that  one  would  expect  to 
measure  with  a  spacecraft  detector.  In  section  3  we  present 
global  MHD,  test  particle,  and  pitch  angle  diffusion  results. 
Finally,  we  conclude  with  a  summary  and  discussion  in 
section  4. 

2.  Global  MHD  Test  Particle  Numerical  Model 
2.1.  Test  Particle  Integrations  in  MHD  Fields 

[7]  An  energetic  particle  distribution  in  the  magneto¬ 
sphere  can  be  modeled  by  following  weighted  Lorentz 
and/or  guiding  center  test  particle  trajectories  in  time- 
dependent  MHD  model  fields;  thus,  the  full  kinetic  effects 
of  particle  interactions  with  MHD  magnetospheric  physics 
is  included.  Since  the  total  energy  density  of  MeV  ions  and 
electrons  in  the  magnetosphere  is  very  small  compared  to 
the  thermal  population,  these  particles  may  be  considered 
noninteracting,  justifying  a  test  particle  approach  [Dessler 
and  Parker ,  1959;  Sckopke,  1966]. 

[s]  Test  particle  trajectories  are  integrated  using  a  fourth- 
order  Runge-Kutta  integrator.  Fast  integration  of  Lorentz 
and  guiding  center  trajectories  is  performed  by  linearly 
interpolating  magnetic  and  electric  fields  from  a  four 
dimensional  (jc,  y ,  z,  /)  grid  to  particle  or  guiding  center 
positions  at  each  time  step.  Since  a  Lorentz  trajectory 
involves  timescales  several  orders  of  magnitude  smaller 
than  a  guiding  center  trajectory,  for  computational  efficiency 
it  is  desirable  to  compute  the  latter.  However,  MeV  particles 
in  the  solar  wind  and  outer  magnetosphere  frequently  do 
not  conserve  their  first  adiabatic  invariant.  In  this  case  it  is 
necessary  to  follow  the  full  Lorentz  motion  of  the  particle, 
and  we  compute  a  particle  trajectory  by  solving  the  rela¬ 
tivistic  Lorentz  equation 


where  v  is  the  velocity  of  a  particle  with  charge  q  and  mass 
m,  7  =  1/^/1  -  v2/c2  is  the  relativistic  factor,  and  c  is  the 
speed  of  light. 

[9]  If  the  gyroperiod  of  a  particle  is  much  smaller  than  the 
timescale  on  which  local  electric  and  magnetic  fields  evolve 
and  the  particle  gyroradius  is  much  smaller  than  the  scale 


3  of  11 


A09215 


KRESS  ET  AL.:  SIMULATIONS  OF  >10  MeV  RADIATION  BELT  ELECTRONS 


A09215 


length  of  local  variations  in  the  fields,  then  the  first 
adiabatic  invariant  is  conserved  by  the  particle  orbit,  and 
we  may  use  the  guiding  center  equations  to  advance  the 
gyroaveraged  center  of  the  particle  position.  The  relativistic 
guiding  center  equations  are 


with  respect  to  a  particle  gyroradius  pgyro,  is  used  to 
determine  the  validity  of  the  guiding  center  approximation: 

|VB| 

c  =  Pgyro^-y-^  (6) 


cb 

Vd  =—  x 

Bq 


—qE  4-  —  Vfl  +  —  (b  V)b 
7  ym  V  / 


7  m 


(2) 


P\\  = 


*Et 


(3) 


Pj 

2  mB 


=  M  =  a  constant. 


(4) 


where  yd  is  the  guiding  center  drift  velocity,  p\\  and  p±  are 
the  parallel  and  perpendicular  components  of  the  relativistic 
particle  momentum,  and  the  first  adiabatic  invariant  M  is  a 
conserved  quantity.  The  brackets  on  the  right-hand  side  of 
equations  (2)  and  (3)  contain  the  electric  field  force  FE,  the 
gradient  drift  force  term  FV£,  the  curvature  force  term  Fcurv, 
and  the  mirror  force  Fmirror  [Northrop,  1963]. 

[10]  When  using  linearly  interpolated  fields,  standard 
adaptive  stepsize  integrators  fail  due  to  lack  of  continuous 
first-  and  higher-order  spatial  derivatives  at  grid  cell  bound¬ 
aries.  For  example,  the  Press  et  al.  [1992]  fourth-order 
Runge-Kutta  quality-controlled  rkqc  routine  uses  the  differ¬ 
ence  between  fourth-  and  fifth-order  approximations  to 
obtain  an  error  approximation  for  adapting  the  time  step. 
An  alternative  method,  used  to  integrate  a  Lorentz  trajectory, 
is  to  set  the  time  step  to  ~  1  / 1 00th  the  instantaneous  particle 
gyroperiod  [e.g.,  Smart  et  al .,  2000].  A  locally  determined 
guiding  center  time  step  is  more  problematic  however,  ( 1 )  due 
to  discontinuities  in  the  linearly  interpolated  magnetic  fields 
at  grid  cell  boundaries,  (2)  since  the  velocity  goes  to  zero  at  a 
mirror  point,  and  (3)  since  the  mirror  force  goes  to  zero  in  the 
magnetic  equatorial  plane.  We  employ  a  guiding  center 
adaptive  stepsize  that  uses  the  instantaneous  guiding  center 
drift  force  terms  to  set  At  at  each  time  step.  A  locally 
determined  adaptive  time  step  for  integrating  the  guiding 
center  equations  is 


At  =  e 


|F£  +  Fvfl  +  Fcurv  +  Fm 


(5) 


where  p  is  the  total  particle  momentum  and  e  is  a  small 
parameter  numerically  determined  to  accurately  integrate  a 
guiding  center  trajectory.  We  find  that  e  ~  0.1  gives 
accurate  results  (e.g.,  £\%  error  for  bounce  and  drift  times 
in  a  dipole  and  to  within  1%  of  a  converged  upon  solution 
in  arbitrary  fields).  In  arbitrary  fields  where  all  the  guiding 
center  force  terms  may  go  to  zero,  e.g.,  in  the  solar  wind,  a 
maximum  time  step  size  should  also  be  used,  e.g..  At  = 
min(A/,  A xmax/Vj). 

[n]  An  adiabaticity  or  “epsilon”  parameter,  which  char¬ 
acterizes  the  length  scale  of  variations  in  the  magnetic  field 


Here  |VB|  is  a  norm  of  the  VB  tensor  (see  Appendix  A). 
This  number  is  used  as  a  switch,  to  go  between  guiding 
center  and  Lorentz  trajectories  when  appropriate.  We  find 
numerically  that  when  £  ^  0.05  there  is  good  agreement 
between  guiding  center  and  Lorentz  methods  over  a  drift 
period.  When  e  £  0.05,  a  trajectory  is  switched  from  a 
guiding  center  to  a  Lorentz  trajectory  using  a  randomly 
chosen  gyrophase  angle.  When  e  ^  0.01,  a  Lorentz 
trajectory  is  switched  to  a  guiding  center  trajectory. 

[12]  In  the  work  presented  here,  all  particles  are  initiated 
with  guiding  center  trajectories,  and  most  remain  in  the 
guiding  center  mode  throughout  the  simulation.  Of  the  few 
that  are  switched  to  the  Lorentz  mode,  most  are  switched 
when  they  encounter  the  magnetopause  and  are  subsequently 
lost  to  the  outer  boundary. 

2.2.  Initializing  the  Test  Particle  Distribution 

[13]  The  computational  domain  for  energetic  particle 
trajectory  tracing  consists  of  a  sphere  centered  on  the  Earth 
with  radius  r  ~  15  RE  (Earth  radii)  located  well  inside  the 
MHD  magnetospheric  model  outer  boundary.  The  dayside 
magnetopause  in  the  MHD  model  is  inside  of  the  sphere, 
while  the  magnetotail  cuts  through  the  sphere  boundary  on 
the  nightside.  Particles  exiting  the  sphere  are  removed  from 
the  simulation.  There  is  also  an  inner  boundary  at  ~2  RE 
which  is  the  inner  boundary  of  the  MHD  magnetospheric 
model  fields.  Particles  striking  this  inner  boundary  are  also 
removed  from  the  simulation. 

[14]  To  initialize  an  outer  radiation  belt  population,  par¬ 
ticles  are  launched  randomly  and  uniformly  from  a  disk  in 
the  equatorial  plane  with  random  uniform  distributions  in 
equatorial  pitch  angle  and  energy.  To  avoid  an  initial  bounce 
phase  bunching,  launch  times  are  distributed  over  a  presi¬ 
mulation  time  interval  greater  than  the  longest  particle 
bounce  period.  Particles  are  launched  with  energies  from 
1  to  7  MeV  and  between  radial  distances  at  3  and  8  RE.  The 
limits  of  the  initial  equatorial  pitch  angle  distribution  are 
determined  by  the  local-time  dependent  loss  cone  in  the 
initial  MHD  model  fields,  and  are  limited  to  £35°  and 
^  145°  by  the  inner  boundary  of  the  MHD  model  at  —2  RE . 

[15]  The  exact  nature  of  the  initial  test  particle  distribu¬ 
tion  is  not  important  except  that  we  wish  to  sample  the 
phase  space  we  are  interested  in  fully  and  in  an  approxi¬ 
mately  uniform  way.  The  initial  distribution  is  subsequently 
weighted,  in  a  postprocessing  step,  using  an  appropriate 
observed  or  model  distribution  which  is  to  be  used  as  the 
initial  condition.  The  particle  weights  are  then  used  to 
obtain  the  fluxes  one  would  expect  to  measure  with  a 
spacecraft  detector.  The  energy  and  radial  limits  of  the 
initial  test  particle  distribution  are  determined  by  the  radi¬ 
ation  belt  model  used  to  weight  the  particles.  In  this  work 
we  use  the  European  Space  Agency  (ESA)  CRRES  radia¬ 
tion  belt  electron  model  [ Vampola ,  1996],  obtained  from 
CRRES  data.  The  model  fluxes  used  are  defined  on  an 
/.-energy  grid  from  L  =  3  to  8  and  E  =  1  to  7  MeV.  Particles 
are  initially  launched  in  this  range.  Any  particles  initially 
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outside  of  these  limits  would  receive  a  particle  weight  of 
zero  and  would  not  contribute  to  the  weighted  particle  flux 
for  the  remainder  of  the  simulation. 

2.3.  Numerically  Determined  Flux 

[16]  Particle  fluxes  are  measured  directly  in  the  code 
using  a  numerical  detector.  The  detector  is  a  disk  in  the 
equatorial  plane.  Test  particle  flux  is  measured  by  counting 
particles  as  they  pass  through  the  disk  and  binning  the 
results  in  position,  energy,  and  equatorial  pitch  angle  space 
(x„  Tj ,  aoAr).  The  flux  integration  time  6t  is  chosen  <7v,min, 
the  smallest  particle  bounce  period,  so  that  we  do  not  over 
sample  a  subset  of  the  trajectories  with  small  bounce 
periods.  For  this  work  fit  ~  0.05  s. 

[17]  The  directional  flux  entering  a  detector  is  in  general  a 
function  of  position  x,  direction  of  incidence  u,  kinetic 
energy  7}  and  time  /,  and  may  be  determined  using 


y(x,u,7\0 


SN 

fiA  cos 0  6Q  fiT  fit' 


(7) 


where  fiN  is  the  number  of  particles  striking  a  surface  of 
area  fiA  with  directions  of  incidence  lying  inside  solid  angle 
fift  oriented  along  the  unit  vector  u,  with  kinetic  energies  in 
the  interval  fiT \  during  the  time  interval  fit.  0  is  the  angle 
between  the  normal  to  fiA  and  the  u  direction  [ Roederer , 
1970,  pp.  85-86].  If  the  distribution  is  uniform  in 
gyrophase,  the  direction  of  incidence  of  the  particles  will 
only  be  a  function  of  pitch  angle  a.  In  this  case  the  solid 
angle  fift  can  be  expressed  as  2n  sin  a  <5a.  If  we  further 
assume  that  in  the  equatorial  plane  a  ^  9,  which  is  usually 
true  to  within  ~1%  for  trapped  particles  in  the  MHD  model 
fields,  the  directional  flux  in  the  ijkth  bin,  with  /,  y,  and  k 
indexing  position,  kinetic  energy,  and  equatorial  pitch  angle, 
respectively,  may  be  obtained  from  the  code  using 


5Z  w" 

Number  ot  llux 

./  -p  \  count,  m  bin  /r»\ 

J  ‘  J  <>l  fiA,  cos  otk  2n  sin  a*  <5a*  <57}  fit  ’ 

where  w„  is  the  particle  weight  (defined  below).  Note  that 
there  is  a  zero  in  the  denominator  for  equatorial ly  mirroring 
particles.  This  difficulty  is  removed  by  choosing  equatorial 
pitch  angle  bins  on  either  side  of  a0  =  90°.  For  simplicity, 
equation  (8)  is  an  expression  for  0th  order  binning  of 
particle  flux  counts,  i.e.,  equivalent  to  the  nearest  grid  point 
(NGP)  method.  A  linear  weighting,  which  also  includes  a 
grid  weighting  factor  for  interpolating  flux  counts  to  (x„  7}, 
<*o)  grid  points,  is  used  in  the  code  to  reduce  noise  in  the 
resulting  distribution  function. 

2.4.  Particle  Weighting 

[is]  As  a  postprocessing  step,  test  particles  are  weighted 
with  an  observed  or  model  flux  distribution  that  is 
chosen  to  serve  as  initial  conditions.  The  test  particle 
flux  counts  are  collected  in  a  file  at  specified  dump 
times  during  the  simulation  and  then  used  along  with 
their  corresponding  particle  weights  to  compute  model 
fluxes  that  would  be  measured  by  a  spacecraft  detector 
throughout  the  simulation. 

[19]  The  particle  weights  are  determined  as  follows.  Near 
the  start  of  the  simulation,  test  particle  fluxes  are  binned  on 


the  same  grid  that  the  weighting  model’s  flux  distribution  is 
specified  on,  to  obtain  an  initial  test  particle  flux  j{csx  parlic|c(/ = 
0).  The  test  particle  flux  is  measured  by  setting  all  particle 
weights  w„  =  1.0  in  equation  (8)  above.  The  initial  test 
particle  flux  is  measured  in  the  code  after  a  time  interval 
greater  than  the  longest  test  particle  drift  period,  in  order  to 
remove  particles  that  were  initiated  in  the  bounce  and  drift 
loss  cones. 

[20]  The  particle  weighting  function  is  determined  by 


w(L  T  n  )  -  7 model (7*/ 1  Tj,aOL) 

yicst  particle  //'»  , 1  —  U ) 


(9) 


When  weighting  particles,  the  weighting  function,  model 
fluxes,  and  initial  test  particle  fluxes  are  all  defined  on  the 
same  grid.  Each  particle  is  given  a  weight  wn  =  w(L„,  7}„ 
c*0 t  =  0).  The  particle  weights  may  be  simply  interpreted 
as  the  number  of  particles  per  test  particle  [e.g.,  see  Hockney 
and  Eastwood ,  1988,  pg.  27].  Each  particle  retains  its 
weight  throughout  the  simulation,  thereby  preserving  the 
total  number  of  particles  along  its  trajectory.  At  later  times 
during  the  simulation,  we  are  not  restricted  to  fluxes  defined 
in  terms  of  the  initial  model  distribution  phase  space 
variables,  e.g.,  particles  with  an  initial  flux  distribution 
defined  in  terms  of  L,  7}  and  a0  may  at  later  times  during  the 
simulation  be  used  to  yield  information  about y‘(x,  7}  a0)  on 
a  different  grid  than  the  one  used  for  weighting. 

[21]  A  number  of  approximations  have  been  made  in  the 
initial  test  particle  weighting  used  in  this  work:  (1)  For 
simplicity,  the  model  L-shell  is  interpreted  as  dipole  L,  i.e. 
radial  distance  in  the  equatorial  plane,  rather  than  using 
Mcllwain  L  determined  in  quiet  Olsen-Pfitzer  model  fields, 
which  was  originally  used  to  sort  CRRES  observations. 
(2)  The  equatorial  pitch  angle  flux  distribution  is  initially 
assumed  to  be  flat  between  35°  and  145°  and  zero  outside 
this  range.  Note  that  one  of  our  present  goals  is  to  isolate  the 
effect  on  the  pitch  angle  distribution  of  particles  accelerated 
by  the  SSC  electric  field  pulse.  Since  outer  zone  electrons 
usually  have  a  maximum  in  their  equatorial  pitch  angle 
distribution  near  90°,  we  expect  the  degree  to  which  the 
final  distribution  is  peaked  near  90°  to  be  underestimated  by 
the  assumption  of  an  initially  flat  pitch  angle  distribution. 
Also,  this  assumption  is  plausible;  e.g.,  see  Seki  et  al.  [2005, 
Figure  5]  which  shows  little  variation  in  ~1  MeV  electron 
fluxes  with  respect  to  equatorial  pitch  angle  near  Z,  ~  4  for 
equatorial  pitch  angles  >45°  and  <135°.  (3)  Most  impor¬ 
tantly,  the  ESA  model  was  produced  using  time  averaged 
CRRES  fluxes.  We  do  not  attempt  to  replicate  the  exact 
conditions  preceding  the  29  October  2003  event  which  are 
largely  unknown  but  rather  explore  radiation  belt  dynamics 
under  extreme  conditions  using  plausible  outer  belt  initial 
conditions. 


3.  Model  Results 

3.1.  MHD  Storm  Sudden  Commencement  Electric 
Field 

[22]  In  the  following  sections  we  present  radiation  belt 
model  results  obtained  by  following  test  particle  trajectories 
in  magnetohydrodynamic  (MHD)  model  fields.  The  mag¬ 
netic  and  electric  fields  are  obtained  from  a  global  MHD 
simulation  of  the  magnetosphere  using  the  Lyon-Feder- 


5  of  11 


A09215 


KRESS  ET  AL.:  SIMULATIONS  OF  >10  MeV  RADIATION  BELT  ELECTRONS 


A09215 


~  30 
co 

g  20 
o 

^  10 

^  -2500 

>  -1500 

f  -500 

>  500 
40.0 

P  13.3 
£ 

£  -13.3 


|  -200 
-500 

0  20  40  60 

Time  (UT  hours  from  0:00  UT  29  Oct  2003) 

Figure  3.  The  29-31  October  2003  ACE  satellite  data 
given  in  SM  coordinates  used  to  drive  the  LFM  magneto- 
spheric  model.  The  Dst  index  from  the  Kyoto  WDC 
Geomagnetic  data  service  is  also  shown.  The  dashed  line  is 
at  the  29  October  2003  SSC.  The  magnetic  and  velocity 
field  components  and  plasma  temperature,  which  are  not 
shown,  are  also  used  as  inputs  to  the  MHD  code. 

Mobarry  (LFM)  code  [Lyon  et  al. ,  2004].  The  solar  wind 
density,  velocity,  temperature  and  magnetic  field  data  from 
the  Anomalous  Composition  Explorer  (ACE)  spacecraft 
were  used  to  drive  the  MHD  model  at  its  sunward  boundary. 
The  29  October  2003  solar  wind  density,  velocity,  and  IMF 
Bz  used  as  inputs  to  the  code  are  shown  in  Figure  3.  A  high- 
resolution  LFM  simulation  was  performed  with  ~325  k  grid 
points  on  a  distorted  spherical  grid.  The  MHD  time  step  was 
~0.02  s.  Field  data  was  dumped  at  ~0.5  s  intervals. 

[23]  An  empirical  plasmaspheric  density  model  [ Gallagher 
et  al .,  1988],  not  usually  included  in  the  MHD  model 
magnetosphere,  was  added  to  the  MHD  fields  immediately 
preceding  the  arrival  of  the  interplanetary  shock  at  ~0600  UT 
on  29  October  2003.  The  LFM  code  does  not  contain  physics 
necessary  to  model  the  evolution  of  the  plasmasphere  on 
longer  timescales  (~hours  to  days),  i.e.,  a  corotation  electric 
field  and  cold  plasma  outflow  source.  The  plasmaspheric 
density  model  is  added  to  modify  the  Alfven  speed  providing 
a  more  accurate  description  of  SSC  pulse  propagation  through 
the  magnetosphere  and  the  subsequent  transient  magneto- 
spheric  oscillations.  The  plasmasphere  does  not  evolve  signif¬ 
icantly  during  the  ~5  s  interval  simulated.  The  propagation  of 
the  SSC  pulse  through  the  inner  magnetosphere  is  sensitive 
to  the  plasmasphere  model  density  profile,  which  is  parame¬ 
terized  in  the  model  solely  with  the  Kp  index.  A  higher  Kp 
reduces  the  plasmaspheric  density  and  moves  the  plasmapause 
earthward  in  the  model.  Kp  is  fixed  for  a  given  300  s  run,  i.e., 
the  plasmasphere  does  not  have  time  to  evolve  significantly  on 
the  SSC  pulse  propagation  timescale.  The  same  plasmasphere 
model  was  used  in  MHD-test  particle  simulations  of  the  March 
1991  event  [Hudson  et  al .,  1997;  Elkington  et  al .,  2002],  with 
resulting  particle  fluxes  in  good  agreement  with  CRRES 
proton  and  electron  measurements. 

[24]  Figure  4  shows  the  azimuthal  component  of  the  SSC 
electric  field  pulse  E0  versus  time  at  several  points  along  the 
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x-axis  in  solar  magnetic  (SM)  coordinates,  i.e.,  approxi¬ 
mately  along  the  day  side  Sun-Earth  line.  The  pulse  is 
shown  for  two  separate  LFM  code  runs  with  two  different 
values  of  the  Kp  input  parameter  into  the  plasmaspheric 
model,  with  Kp-  3  (Figure  4a)  and  with  Kp  =  4  (Figure  4b). 
In  each  case,  the  initial  E0  displacement  is  in  the  westward 
(negative  0)  direction,  which  can  be  understood  by 
considering  Faraday’s  law  and  a  compression  of  the  Earth’s 
northwardly  directed  magnetic  field.  Thus  the  leading 
portion  of  the  bipolar  pulse  shown  in  Figure  4  causes 
electrons  and  ions  to  E  x  B  drift  inward  in  L-shell.  Particles 
in  drift  resonance  with  the  azimuthal  propagation  of  the 
pulse  may  be  transported  earthward  several  Earth  radii, 
undergoing  significant  adiabatic  energization.  A  rough 
estimate  for  the  azimuthal  propagation  velocity  of  the  pulse 
can  be  obtained  by  considering  the  time  between  the 
appearance  of  the  pulse  on  the  dayside  at  +3  RE  and  on  the 
nightside  at  -3  RE.  In  the  case  with  Kp  =  3  this  yields  37 tRe/ 
60s  «  1000  km/s.  In  general,  the  effect  of  the  addition  of  the 
plasmaspheric  density  is  to  slow,  steepen,  and  intensify  the 
SSC  electric  field  pulse  in  the  inner  magnetosphere.  On  the 
day  preceding  the  arrival  of  the  29  October  2003  CME,  the 
Kp  index  fluctuated  between  Kp  =  3  and  Kp  =  5  (Kyoto 
WDC  Geomagnetic  data  service).  On  the  morning  of 
29  October  2003,  preceding  the  arrival  of  the  shock,  Kyoto 
WDC  reported  Kp  =  4.  Test  particle  results  in  fields  from 
MHD  runs  using  Kp  =  3  and  Kp  -  4  are  presented  below. 
The  case  with  Kp  =  3  results  in  higher  flux  levels  of 
>10  MeV  electrons  inside  of  L  ~  3  than  the  Kp  =  4  case. 
MHD  simulations  with  Kp  =  2  and  with  no  plasmaspheric 
density  model  were  also  performed  (not  shown).  These 
additional  runs  follow  the  trend  illustrated  in  Figure  4,  i.e., 
an  increased  density  in  the  inner  magnetosphere  produces  a 
slower,  narrower,  and  larger  SSC  electric  field  pulse. 
The  case  with  no  plasmasphere  model  included  in  the 
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Figure  4.  Azimuthal  component  of  the  MHD  model  SSC 
electric  field  pulse  versus  time  at  several  points  along  the  x- 
axis  in  solar  magnetic  (SM)  coordinates,  i.e.,  approximately 
along  the  Earth-Sun  line.  Results  from  two  separate  MHD 
magnetospheric  model  code  runs  are  shown,  with  (a)  Kp  =  3 
and  (b)  Kp  -  A  used  as  input  parameters  into  the 
plasmaspheric  density  model. 
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Figure  5.  Time  snapshot  of  the  azimuthal  component  of 
the  MHD  model  SSC  electric  field  pulse  in  the  equatorial 
plane.  The  dashed  line  shows  the  trajectory  of  a  single 
adiabatically  accelerated  guiding  center  electron  in  drift 
resonance  with  the  pulse  as  it  propagates  from  the  dayside 
to  nightside.  The  initial  and  final  energies  of  the  particle  are 
~5  and  15  MeV,  respectively. 

MHD  density  produces  a  maximum  electric  field  pulse  of 
~60  mV/m  in  the  inner  magnetosphere,  which  is  not 
sufficient  to  produce  a  significant  >10  MeV  belt.  The  run 
with  Kp  =  2  produces  a  maximum  electric  field  amplitude 
~  1 50  mV/m,  only  slightly  larger  than  the  run  with  Kp  =  3. 

[25]  A  time  snapshot  of  E0  in  the  equatorial  plane  is 
shown  in  Figure  5.  Also  shown  in  Figure  5  is  the  trajectory 
of  a  single  guiding  center  electron  that  is  in  drift  resonance 
with  the  pulse,  moving  with  the  crest  of  the  pulse  as  it 
propagates  from  the  dayside  to  nightside.  The  initial  and 
final  energies  of  the  particle  are  ^5  and  15  MeV, 
respectively.  The  electron  trajectory  shown  is  equatorial ly 
mirroring  (equatorial  pitch  angle  is  90°).  A  meridional  plot 
of  the  electric  field  magnitude  (not  shown)  reveals  that  the 
SSC  electric  field  pulse  is  mainly  in  the  equatorial  plane, 
thus  preferentially  accelerating  equatorial  ly  mirroring  par¬ 
ticles  which  remain  in  the  strongest  portion  of  the  pulse. 
This  effect  is  enhanced  by  an  additional  focusing  of  the 
pulse  into  the  equatorial  plane  as  it  enters  the  inner 
magnetosphere. This  result  is  illustrated  in  Figure  6  which 
shows  the  SSC  E0  at  two  separate  time  snapshots  as  it  enters 
the  inner  magnetosphere  plotted  in  the  noon-midnight 
meridional  plane  in  solar  magnetic  (SM)  coordinates.  The 
pulse  is  focused  into  the  equatorial  plane  as  it  enters  the 
inner  magnetosphere.  In  the  inner  magnetosphere,  the  pulse 
remains  near  the  equatorial  plane  as  it  propagates  toward  the 
nightside. 

3.2.  Test  Particle  Model  Results 

[26]  In  each  run,  ~2.4  million  test  particle  trajectories  are 
computed  during  an  ~5  min  interval  from  the  MHD 
simulation  that  includes  the  initial  impact  of  the  interplan¬ 
etary  shock  on  the  magnetosphere  at  ~0600  UT,  propaga¬ 
tion  of  the  resulting  fast  mode  magnetosonic  pulse  through 


the  magnetosphere,  and  during  several  large  transient  ULF 
oscillations  following  the  arrival  of  the  shock  shown  in 
Figure  4.  An  overview  of  the  weighted  test  particle  radiation 
belt  results  from  the  Kp  =  3  case  are  shown  in  Figure  7. 
Omnidirectional  integrated  >10  MeV  fluxes  are  plotted  in 
the  equatorial  plane  in  SM  coordinates  at  four  snapshots  in 
time  from  the  simulation.  Since  the  initial  model  distribu¬ 
tion  has  no  particles  with  energies  above  7  MeV,  there  are 
initially  zero  >10  MeV  fluxes.  The  initially  localized 
injection  appears  at  ~1500  local  time.  Figure  7  nicely 
illustrates  the  source  of  the  drift  echoes  observed  by  a 
spacecraft  particle  detector  [e.g.,  Blake  et  al.,  1992, 
Figure  1],  i.e.,  a  sudden  appearance  of  heightened  fluxes, 
with  higher  energy  particles  reaching  the  detector  before 
lower  energies  due  to  a  S7B  drift  velocity  dispersion  of  the 
initially  localized  injection.  At  a  fixed  location  there  is  a 
gradual  rise  in  fluxes  as  lower  energies  reach  the  detector 
until  a  sudden  drop  in  fluxes  occurs  when  the  detector's 
lower-energy  cutoff  is  reached.  Subsequent  drift  echoes  are 
gradually  diminished  by  energy  dispersion. 

[27]  Figure  8  shows  10  MeV  equatorial  pitch  angle 
distributions  in  the  newly  formed  belt  ~5  min  after  storm 
onset.  To  produce  the  distributions,  weighted  particle  fluxes 
are  binned  in  dipole  /.-shell  (radial  distance  in  the  equatorial 
plane),  energy,  and  equatorial  pitch  angle.  In  the  case  with 
Kp  =  3  the  peak  in  10  MeV  flux  occurs  at  L  ~  3.0,  with  the 
corresponding  distribution  shown  in  Figure  8b.  Figure  8a 
shows  the  10  MeV  equatorial  pitch  angle  distribution  inside 
the  flux  peak,  at  Z,  ~  2.5.  Figures  8c  and  8d  show 
distributions  resulting  from  the  run  with  Kp  =  4,  at  the 


29  Oct  2003  MHD  SSC  (mV/m) 


Figure  6.  MHD  EQ  peak  at  6  RE  and  at  3  RE  (two  separate 
time  snapshots)  on  the  dayside  in  the  noon-midnight 
meridional  plane  in  solar  magnetic  (SM)  coordinates. 
Contours  from  —80  to  —140  mV/m  in  steps  of  -20  mV/ 
m  are  shown.  Dashed  lines  show  magnetic  field  lines 
through  x  =  6  Re  and  jc  =  3  RE  traced  in  the  MHD  fields  in 
each  respective  time  snapshot. 
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Figure  7.  Integral  omnidirectional  fluxes  of  >10  MeV  electrons  in  the  equatorial  plane  in  solar  magnetic 
(SM)  coordinates  at  four  snapshots  in  time  from  the  MHD  test  particle  simulation,  from  the  MHD 
magnetospheric  model  run  using  Kp  =  3  as  an  input  parameter  to  the  plasmaspheric  density  model.  The 
initial  test  particle  distribution  is  weighted  using  the  European  Space  Agency  (ESA)  CRRES  radiation 
belt  electron  model  [Vampola,  1996].  The  model  fluxes  used  for  initial  conditions  are  defined  on  an 
L-energy  grid  from  L  =  3  to  8  and  E  =  1  to  7  MeV;  thus  the  initial  >10  MeV  flux  in  the  model  is  zero. 


location  of  the  10  MeV  peak  at  L  ~  3.5,  and  inside  the 
10  MeV  peak  at  L  ~  3.  In  each  case,  the  solid  line  plot  is  a 
least  squares  fit  to  A  sin"(a„).  The  resulting  distributions  are 
strongly  peaked  around  90°  with  fluxes  falling  to  zero 
near  ±60°.  Also  note  that  the  distributions  become  more 
peaked  with  decreasing  L.  The  line  plot  shown  in  Figure  8b 
with  Kp  =  3  at  3.0  RE  is  used  as  the  initial  condition  for  the 
pitch  angle  diffusion  calculation  presented  in  the  next  section. 

3.3.  Pitch  Angle  Diffusion 

[28]  An  estimate  for  the  timescale  of  the  pitch  angle 
diffusion  of  an  equatorially  mirroring  population  of  radia¬ 
tion  belt  electrons,  to  a  spacecraft  position  in  low  Earth 
orbit,  may  be  obtained  by  solving  a  one-dimensional  pitch 
angle  diffusion  equation.  Bounce-averaged  equatorial  pitch 
angle  diffusion  coefficients  are  calculated  according  to 
quasi-linear  theory  in  the  high-density  approximation 
[Lyons,  1974a,  1974b],  using  the  computational  techniques 
of  Albert  [1999]  and  the  parameters  of  Meredith  et  al. 
[2006].  These  values  are  based  on  CRRES  data  and  suc¬ 
cessfully  reproduced  the  observed  decay  times  of  1.09  MeV 
electrons  at  3  <  L  <  4  and  over  a  broader  range  in  L  for 
lower  energies. 

[29]  The  model  has  a  density  ratio  fpelfce  =  8.9 
(corresponding  to  ne  =  1086.6)  and  a  hiss  amplitude  of 
34.5  pT.  The  power  spectral  density  is  a  truncated  Gaussian 


peaked  at  550  Hz  with  width  300  Hz  and  lower  and  upper 
cutoffs  at  100  Hz  and  2000  Hz,  respectively.  The  “small 
wavenormal  model”  is  also  a  Gaussian  (in  x  -  tan#),  peaked 
at  x  =  0  with  width  xw  =  tan  20°,  truncated  at  ±rmax  =  tan 
30°.  Since,  for  high-energy  particles,  cyclotron  resonance 
can  occur  with  large  harmonic  numbers,  n  up  to  ±100  was 
kept,  although  most  large  n  values  are  quickly  eliminated  by 
the  computational  procedure. 

[30]  The  solid  curve  in  Figure  9  shows  the  resulting 
diffusion  coefficients  for  10  MeV  electrons,  while  the 
dashed  curve  shows  1  MeV  values  for  comparison.  The 
10  MeV  values  were  used  to  evolve  the  electrons  in  time 
according  to  the  one-dimensional  pitch  angle  diffusion 
equation 

-77-  =  ~  : - ! - ( D0olXs  T  sin  ot()  cos  a„  Y  (10) 

Ot  T  sin  aQ  cos  a0  da„  \  0ao) 

where  T  is  the  bounce  period,  with  boundary  conditions 
f(aOL()  =  0  and  ( djJdaa\Xo  =  90°  =  0.  The  diffusion 
simulation  was  performed  at  L  =  3  in  a  dipole,  where  the 
equatorial  pitch  angle  at  the  edge  of  the  loss  cone  is  at  ~8°. 
The  initial  distribution  was  modeled  as  j  =  A  sin"  aQ  with  A  = 
40,  n  =  22  (Figure  8b).  Grid  resolution  was  1  degree,  with 
timestep  10  4  days.  The  curves  shown  in  Figure  10  are 
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Figure  8.  Equatorial  pitch  angle  distributions  of  10  MeV 
differential  flux  for  various  Kp  and  radial  distances.  The 
asterisks  show  nonzero  fluxes  obtained  with  equation  (8) 
with  28  a0<  bins  uniformly  spaced  in  cos  a0  between  0  and 
180  deg.  The  solid  line  is  a  least  squares  fit  to  A  sin"(a0) 
yielding:  (a)  A  =  3.2,  n  =  40;  (b)  A  =  40,  n  =  22;  (c)  A  =  12, 
n  =  30;  (d)  A  =  20,  n  =  18.  Note  that  in  each  figure,  the  axis 
has  been  scaled  to  span  2  decades  and  the  maximum  value 
of  the  flux  has  been  located  the  same  distance  below  the  top 
of  the  vertical  axis  so  that  the  rate  the  fluxes  fall  off  from 
their  maximum  values  can  be  compared. 


snapshots  over  200  days  at  increments  of  10  days,  starting 
from  the  dark  blue  curve.  The  flux  takes  1  to  2  months  to 
be  discernable  at  the  low  equatorial  pitch  angles 
observable  by  SAMPEX  at  ~10  degrees  (dashed  line) 
and  3  to  4  months  for  the  steady  state  pitch  angle  profile 
to  be  established,  which  is  in  agreement  with  the  timescale 
for  the  delay  in  the  appearance  of  peak  flux  levels  at 
SAMPEX  seen  in  Figure  1. 

4.  Summary  and  Discussion 

[31]  At  energies  in  the  10s  of  MeV  range  the  structure  of 
the  inner  zone  radiation  belts  is  largely  shaped  by  a  few 
geomagnetic  storms  driven  by  high-speed  interplanetary 
shock  compressions  of  the  magnetopause.  The  MHD  test 
particle  model  results  show  that  the  29  October  2003  SSC 
electric  field  pulse  produces  a  new  belt  of  >10  MeV 
electrons  inside  of  L  ~  3  with  an  average  quiet-time  outer 
belt  model  source  population  assumed.  The  newly  formed 
>10  MeV  electron  belt  has  its  equatorial  pitch  angle  distri¬ 
bution  strongly  peaked  near  90°  and  becomes  more  peaked 
with  lower  L  as  shown  in  Figure  8.  There  are  two  primary 
reasons  for  the  resulting  peaked  equatorial  pitch  angle 
distribution  in  the  new  belt:  (1)  The  electrons  are  acceler¬ 
ated  through  conservation  of  the  1st  adiabatic  invariant 
perpendicular  to  the  magnetic  field  increasing  pjp\\  and 
bringing  their  equatorial  pitch  angles  closer  to  90°,  and 
(2)  the  SSC  electric  field  pulse  is  predominantly  in  the 
equatorial  plane,  preferentially  accelerating  equatorially 
mirroring  particles  that  spend  more  time  in  the  pulse.  The 


parameters  of  Horne  et  al.  [2005a,  2005b]  (See  text  for 
model  parameters).  The  solid  curve  shows  the  resulting 
diffusion  coefficients  for  10  MeV  electrons,  while  the 
dashed  curve  shows  1  MeV  values  for  comparison. 

SSC  electric  field  pulse  is  focused  into  the  equatorial  plane 
as  it  propagates  into  the  plasmasphere,  which  has  been 
added  to  the  MHD  fields  immediately  before  the  arrival  of 
the  shock  using  a  /^-dependent  empirical  density  model.  In 
general  the  effect  of  including  a  plasmaspheric  density 
model  in  the  MHD  simulations  is  to  enhance  the  amplitude 
of  the  SSC  electric  field  pulse  and  decrease  its  speed  in  the 
inner  magnetosphere.  It  is  necessary  to  include  a  realistic 
plasmaspheric  density  in  the  MHD  magnetospheric  model 
to  produce  an  SSC  electric  field  pulse  large  enough  to 
transport  electrons  over  several  L-shells  producing  a 
significant  >10  MeV  belt.  This  is  consistent  with  Gannon 
et  al.  [2005]  who  find  that  an  SSC  E-field  pulse  ~10  mV/m 


Figure  10.  Equatorial  pitch  angle  distribution  of  10  MeV 
electrons  evolved  using  (10)  over  200  days  at  increments  of 
10  days,  starting  from  the  dark  blue  curve.  The  flux  takes 
1  to  2  months  to  be  discernable  at  the  low  equatorial  pitch 
angles  observable  by  SAMPEX  at  ~10  degrees  (dashed 
line)  and  3  to  4  months  for  the  steady  state  pitch  angle 
profile  to  be  established. 
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is  more  than  an  order  of  magnitude  too  small  to  produce  a 
significant  new  belt  with  energies  ~10  MeV.  The  delay, 
between  the  injection  of  >10  MeV  electrons  on  29  October 
2003  and  the  appearance  of  peak  flux  levels  at  SAMPEX 
seen  in  Figure  1 ,  indicates  a  pitch  angle  diffusion  timescale  ~ 
months  which  is  in  agreement  with  an  estimate  obtained 
by  solving  a  one-dimensional  pitch  angle  diffusion 
equation. 

[32]  On  the  basis  of  our  results,  we  expect  the  location, 
energy,  and  equatorial  pitch  angle  distribution  of  the  newly 
formed  belt  to  be  dependent  on  the  state  of  the  outer  zone 
electrons  immediately  preceding  the  SSC.  In  this  work,  we 
have  not  attempted  to  reproduce  the  exact  state  of  the 
plasmasphere  or  outer  electron  belt  preceding  the  29  October 
2003  event.  Our  main  goal  has  been  to  study  radiation  belt 
dynamics  under  extreme  conditions  using  a  plausible  initial 
state  for  the  outer  belt. 

[33]  It  is  not  possible  to  make  a  direct  comparisons  with 
observation  in  this  study,  of  the  location  and  intensity  of  the 
newly  formed  radiation  belt,  since  there  was  no  satellite  in 
near  equatorial  orbit  capable  of  sampling  the  ultrarelativistic 
electrons  as  CRRES  did  during  the  24  March  1991  event. 
However,  the  model  does  not  produce  a  peak  in  >10  MeV 
electron  flux  near  L  =  2,  where  the  new  belt  appears  in  the 
SAMPEX  summary  plot  (Figure  1)  several  months  follow¬ 
ing  the  29  October  2003  storm.  For  example,  in  the  case 
with  Kp  —  3  used  as  an  input  parameter  to  the  plasmaspheric 
density  model,  the  peak  in  10  MeV  flux  forms  at  ~3  RE. 
Peak  fluxes  at  higher  energies  occur  at  lower  L  until  the 
inner  boundary  of  the  MHD  magnetospheric  model  is 
reached  at  ~2  RE.  The  electric  field  pulse  is  modified  in 
this  region  by  the  conditions  imposed  at  the  inner  boundary 
on  the  MHD  fields:  the  normal  component  of  the  magnetic 
flux  through  the  inner  boundary  surface  is  held  constant, 
and  the  normal  component  of  the  velocity  at  the  inner 
boundary  is  zero,  greatly  diminishing  the  inductive  electric 
field  near  the  inner  boundary  and  thereby  reducing  inward 
transport  due  to  E  x  B  drift  near  L  ~  2.  The  Courant 
condition  on  pulse  propagation  across  a  grid  cell  inhibits 
moving  the  inner  boundary  closer  to  Earth,  where  the 
AlfVen  speed  is  greater  [Lyon  et  al ,  2004]. 

[34]  While  progress  has  been  made  with  2-D  guiding 
center  particle  simulations  modeling  the  24  March  1991 
injection,  for  which  near  equatorial  plane  measurements 
were  available,  we  have  since  relied  mainly  on  continuous 
low  altitude  polar  orbiting  spacecraft  measurements  in  and 
near  the  loss  cone.  A  number  of  storms  subsequent  to  the 
March  1991  event  show  evidence  of  prompt  shock-drift 
injections.  Blake  et  al  [2005]  examined  the  geoeffective¬ 
ness  of  shocks  populating  the  radiation  belts  and  found  a 
correlation  with  sharp  increase  in  the  H-component,  ground 
magnetometer  perturbation  which  is  the  signature  of  an  SSC 
event.  The  March  1991  event,  two  in  November  2001 
which  produced  trapping  of  solar  energetic  protons  [Kress 
et  al .,  2004,  2005;  Hudson  et  al .,  2004;  Mazur  et  al .,  2006], 
and  both  the  21  February  1994  and  28  October  2003  events 
appear  to  meet  the  criterion  of  producing  a  large  inductive 
electric  field  pulse  capable  of  transporting  earthward  to  L  < 

3  particles  with  azimuthal  drift  velocities  comparable  to  the 
pulse  propagation  speed  (~1000  km/s).  CRRES  electron 
data  following  the  24  March  1991  storm  shows  evidence  of 
pitch  angle  distribution  strongly  peaked  around  90°  in  the 


newly  formed  >10  MeV  belt.  A  delay  ~  months  before  the 
appearance  of  peak  fluxes  at  LEO  following  the  2 1  February 
1994  and  29  October  2003  storms  suggest  similarly  peaked 
distributions,  illustrating  the  need  for  an  understanding  of 
pitch  angle  evolution  to  model  energetic  particle  fluxes  in  the 
radiation  belts. 

[35]  The  model  results  show  that  a  quiet  time  outer  zone 
electron  distribution  provides  a  source  population  for  the 
observed  injection,  however,  high  fluxes  of  solar  energetic 
particles  were  also  present  at  the  arrival  of  the  interplanetary 
shock  on  29  October  2003,  with  flux  levels  of  solar 
energetic  electrons  in  the  1-7  MeV  range  approaching 
those  usually  observed  in  the  outer  belts  [Mewaldt  et  al ., 
2005].  Solar  energetic  electrons  (SEEs)  provide  a  second 
possible  source  population  for  the  newly  formed  10 
20  MeV  electron  belt.  SEEs  with  a  relativistic  gamma  factor 
£10  have  gyroradii  approaching  the  scale  lengths  of  mag¬ 
netic  field  gradients  in  the  outer  radiation  belts,  suggesting 
they  may  be  promptly  trapped  in  the  magnetosphere  as  has 
previously  been  shown  for  solar  energetic  ions  [Kress  et  al ., 
2005;  Mazur  et  al .,  2006].  We  reserve  for  future  investiga¬ 
tion  the  role  of  SEEs  as  a  source  population  for  trapped 
>10  MeV  electrons  in  the  inner  zone. 


Appendix  A:  Nonadiabaticity 

[36]  The  adiabaticity  or  “epsilon”  parameter,  used  to 
determine  the  validity  of  the  guiding  center  approximation, 
is  frequently  expressed 


e  =  max(Jpgyro  — 


(Al) 


However,  in  arbitrary  three-dimensional  fields  this  expres¬ 
sion  does  not  take  into  account  magnetic  rotational  shear. 
As  a  simple  example,  consider  the  magnetic  field 


B  =  sin(rtz)Jc  +  cos(/ir)v,  (A2) 

for  which  e  =  0  for  arbitrarily  rapid  variation  in  the  z 
direction  (i.e.,  arbitrarily  large  n).  The  field  expressed  by 
equation  (A2)  has  no  curvature  or  gradient  in  its  magnitude, 
however,  there  is  a  rotational  shear  as  we  proceed  in  the  z 
direction.  Rotational  shear  is  not  uncommon  in  the 
magnetosphere  and  is  produced  at  any  place  we  find  a 
current  sheet;  therefore  in  general  we  suggest  a  norm  of  the 
VB  tensor: 


which  captures  all  components  of  variations  in  the  magnetic 
field  (D.  C.  Montgomery,  private  communication,  2006). 
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